Excess mortality: combining three countries
Data
expected_deaths_monthly <- read_rds("data/expected_deaths_monthly.Rds") %>%
mutate(Pandemic = case_when(
(Year >= 1889 & Year <= 1892) ~ 1890,
(Year >= 1917 & Year <= 1920) ~ 1918,
(Year >= 1956 & Year <= 1959) ~ 1957,
(Year >= 2019 & Year <= 2020) ~ 2020
))Figure 1 - diffs
Monthly
Yearly
Figure 2 - four pandemics
Absolute
focus_years <- c(1890, 1918, 1957, 2020)
# plot_years <- c(seq(1890 - 1, 1890 + 2),
# seq(1918 - 1, 1918 + 2),
# seq(1957 - 1, 1957 + 2),
# seq(2020 - 1, 2020 + 2))
for(country in unique(expected_deaths_monthly$Country)) {
print(paste0("Plotting ", country))
max <- expected_deaths_monthly %>%
filter(Country == country) %>%
filter(!is.na(Pandemic)) %>%
summarise(deaths = max(deaths),
lpi_flu_hc = max(lpi_flu_hc),
upi_flu_hc = max(upi_flu_hc)) %>%
rowwise() %>%
mutate(max = max(c(deaths, lpi_flu_hc, upi_flu_hc)))
print(paste0("Max for x is ", max$max))
for(year in focus_years) {
print(paste0(" Year ", year))
plot_data <- expected_deaths_monthly %>%
filter(Country == country) %>%
filter(Pandemic == year)
name <- paste(country, year, sep = "_")
if (year == 2020) {
my_title = paste0(country, ": ",
as.character(year - 1), "-", as.character(year))
} else {
my_title = paste0(country, ": ",
as.character(year - 1), "-", as.character(year + 2))
}
plot <- plot_data %>%
ggplot(aes(x = Date)) +
geom_rect(xmin = as.Date(paste0(year, "-01-01")),
xmax = as.Date(paste0(year, "-12-31")),
ymin = -Inf, ymax = Inf,
fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
geom_ribbon(aes(ymin = lpi_flu_hc, ymax = upi_flu_hc), alpha = 0.1) +
geom_line(aes(y = fit_flu_hc), colour = "grey50") +
# geom_line(aes(y = median_flu_year), colour = "#2c7bb6") +
geom_line(aes(y = deaths), color = "#d7191c", size = 0.5) +
scale_x_date(limits = c(min(plot_data$Date), max(plot_data$Date)),
date_breaks = "6 months", date_labels = "%b") +
xlab("") +
# ylab("Observed, 10y median & predicted deaths") +
ylab("Observed & predicted deaths") +
# labs(caption = "Excludes extreme years") +
ggtitle(my_title) +
scale_y_continuous(limits = c(0, max$max),
labels = label_number(accuracy = 0.1, suffix = "K",
scale = 1e-4)) +
theme_bw() +
theme(plot.title = element_text(size = 10),
axis.text.x = element_text(size = 8),
axis.title.y = element_text(size = 9),
axis.text.y = element_text(size = 8))
assign(name, plot)
}
}[1] "Plotting Switzerland"
[1] "Max for x is 10766"
[1] " Year 1890"
[1] " Year 1918"
[1] " Year 1957"
[1] " Year 2020"
[1] "Plotting Sweden"
[1] "Max for x is 17278"
[1] " Year 1890"
[1] " Year 1918"
[1] " Year 1957"
[1] " Year 2020"
[1] "Plotting Spain"
[1] "Max for x is 163422"
[1] " Year 1890"
[1] " Year 1918"
[1] " Year 1957"
[1] " Year 2020"
plot_spacer() + Spain_1918 + Spain_1957 + Spain_2020 +
Sweden_1890 + Sweden_1918 + Sweden_1957 + Sweden_2020 +
Switzerland_1890 + Switzerland_1918 + Switzerland_1957 + Switzerland_2020 +
plot_layout(widths = rep(c(2, 2, 2, 1), 3)) +
plot_layout(ncol = 4) ggsave("paper/four_pndemics.png", dpi = 300, width = 297, height = 210, units = "mm")Percent
plot_years <- c(seq(1890 - 1, 1890 + 2),
seq(1918 - 1, 1918 + 2),
seq(1957 - 1, 1957 + 2),
seq(2020 - 1, 2020 + 2))
plot_data <- expected_deaths_monthly %>%
filter(Year %in% plot_years) %>%
mutate(Percent = ((deaths - fit_flu_hc) / fit_flu_hc),
Difference = ifelse(deaths > fit_flu_hc, "More", "Less")) %>%
group_by(Country, Pandemic) %>%
mutate(time = row_number()) %>%
ungroup()
plot_data %>%
ggplot(aes(x = time)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = 13, colour = "grey80") +
geom_vline(xintercept = 24, colour = "grey80") +
geom_rect(xmin = as.Date(paste0(year, "-01-01")),
xmax = as.Date(paste0(year, "-12-31")),
ymin = -Inf, ymax = Inf,
fill = "#ffffbf", alpha = 0.01, inherit.aes = FALSE) +
geom_col(aes(y = Percent, fill = Difference), size = 0.5) +
facet_grid(vars(Country), vars(Pandemic), scales = "free_x") +
scale_x_continuous(limits = c(min(plot_data$time)-1, max(plot_data$time)+1),
breaks = c(1, 13, 24, 36, 48),
labels = c("-12m", "Start", "End", "+12m", "+24m")) +
xlab("") +
ylab("% above predicted deaths") +
ggtitle("Excess deaths as % of observed") +
scale_y_continuous(labels = percent) +
scale_fill_manual(values = c("#67A9CF", "#EF8A62")) +
guides(fill = guide_legend(reverse = TRUE)) +
theme_bw() +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(),
panel.grid.minor.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "grey40")) ggsave("paper/four_pndemics_perc.png", dpi = 300, width = 297, height = 210, units = "mm")Table 1
pandemic_year <- c("1890","1918","1957","2020")
table_1 <- expected_deaths_monthly %>%
select(Country,Year, deaths, pop, fit_flu_hc) %>%
filter(Year %in% pandemic_year)%>%
group_by(Year, Country) %>%
summarise(Expected_Mortality = round(sum(fit_flu_hc), 0),
Deaths_num = sum(deaths),
Cum_excess_death = round(sum(deaths)-sum(fit_flu_hc), 0),
Percent_excess_death = round((Cum_excess_death/Expected_Mortality)*100, 1))%>%
ungroup()
openxlsx::write.xlsx(table_1,"paper/table1.xlsx", row.names = FALSE)| Year | Country | Expected_Mortality | Deaths_num | Cum_excess_death | Percent_excess_death |
|---|---|---|---|---|---|
| 1890 | Sweden | 73077 | 81824 | 8747 | 12.0 |
| 1890 | Switzerland | 58397 | 61805 | 3408 | 5.8 |
| 1918 | Spain | 451779 | 695758 | 243979 | 54.0 |
| 1918 | Sweden | 78642 | 104591 | 25949 | 33.0 |
| 1918 | Switzerland | 50473 | 75034 | 24561 | 48.7 |
| 1957 | Spain | 282888 | 293502 | 10614 | 3.8 |
| 1957 | Sweden | 70064 | 73132 | 3068 | 4.4 |
| 1957 | Switzerland | 52280 | 51066 | -1214 | -2.3 |
| 2020 | Spain | 425484 | 494604 | 69120 | 16.2 |
| 2020 | Sweden | 88975 | 95411 | 6436 | 7.2 |
| 2020 | Switzerland | 67893 | 75971 | 8078 | 11.9 |
Appendix
## Crude
expected_deaths_monthly %>%
ggplot(aes(x = Date, y = deaths)) +
geom_col() +
geom_rect(xmin = as.Date(paste0(1890, "-01-01")),
xmax = as.Date(paste0(1890, "-12-31")),
ymin = 0, ymax = Inf,
fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
geom_rect(xmin = as.Date(paste0(1918, "-01-01")),
xmax = as.Date(paste0(1918, "-12-31")),
ymin = 0, ymax = Inf,
fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
geom_rect(xmin = as.Date(paste0(1957, "-01-01")),
xmax = as.Date(paste0(1957, "-12-31")),
ymin = 0, ymax = Inf,
fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
geom_rect(xmin = as.Date(paste0(2020, "-01-01")),
xmax = as.Date(paste0(2020, "-12-31")),
ymin = 0, ymax = Inf,
fill = "#d7191c", alpha = 0.01, inherit.aes = FALSE) +
facet_grid(vars(Country), scales = "free_y") +
scale_y_continuous(labels = label_number(accuracy = 1L, suffix = "K", scale = 1e-4)) +
theme_minimal() +
xlab("") + ylab("Deaths") +
ggtitle("Monthly number of deaths")Per population
expected_deaths_monthly %>%
ggplot(aes(x = Date, y = (deaths/pop)*10000, colour = Country)) +
geom_line() +
theme_minimal() +
xlab("") + ylab("Deaths / 10,000 population") +
ggtitle("Monthly number of deaths")